require(stats)
require(graphics)
#library("ggplot2")
require(ggplot2)
#library("grDevices")
require(grDevices)
#library("gridExtra")
require(gridExtra)
#library("grid")
require(grid)
#install.packages("Hmisc")
require(Hmisc)

# Please make sure to change the path to where the appropriate dataset is located on your computer

# FR
fr.k1.data <- read.csv("PNAS_VARILL_FR_K1.csv")
fr.k2.data <- read.csv("PNAS_VARILL_FR_K2.csv")
fr.k3.data <- read.csv("PNAS_VARILL_FR_K3.csv")
# ST
st.k1.data <- read.csv("PNAS_VARILL_ST_K1.csv")
st.k2.data <- read.csv("PNAS_VARILL_ST_K2.csv")
st.k3.data <- read.csv("PNAS_VARILL_ST_K3.csv")
# FAM
fam.k1.data <- read.csv("PNAS_VARILL_FAM_K1.csv")
fam.k2.data <- read.csv("PNAS_VARILL_FAM_K2.csv")
fam.k3.data <- read.csv("PNAS_VARILL_FAM_K3.csv")
#LS
ls.k1.data <- read.csv("PNAS_VARILL_LS_K1.csv")
ls.k2.data <- read.csv("PNAS_VARILL_LS_K2.csv")
ls.k3.data <- read.csv("PNAS_VARILL_LS_K3.csv")
#HS
hs.k1.data <- read.csv("PNAS_VARILL_HS_K1.csv")
hs.k2.data <- read.csv("PNAS_VARILL_HS_K2.csv")
hs.k3.data <- read.csv("PNAS_VARILL_HS_K3.csv")
#O
open.k1.data <- read.csv("PNAS_VARILL_O_K1.csv")
open.k2.data <- read.csv("PNAS_VARILL_O_K2.csv")
open.k3.data <- read.csv("PNAS_VARILL_O_K3.csv")
#NO
no.k1.data <- read.csv("PNAS_VARILL_NO_K1.csv")
no.k2.data <- read.csv("PNAS_VARILL_NO_K2.csv")
no.k3.data <- read.csv("PNAS_VARILL_NO_K3.csv")
#
df.list <- list(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data,
                no.k2.data,fr.k2.data,fam.k2.data,ls.k2.data,open.k2.data,hs.k2.data,st.k2.data,
                no.k3.data,fr.k3.data,fam.k3.data,ls.k3.data,open.k3.data,hs.k3.data,st.k3.data)

sponsor.require <- 1
high.skilled.p <-  2
low.skilled.p <-  3
p.tourist.visa <- 4
p.fakedocs <-  5
no.econ.stud <- 6
p.all.illegal <- 7
k <- 8
my.seed <- 9

legaltourist.attempts <- c(10:29)

ability.legal.cols <- c(30:90)
ability.fakedocs.cols <- c(91:151)
ability.work.cols <- c(152:212)

legal.cols <- c(274:334)
fi.cols <- c(335:395)
sl.cols <- c(396:456)
stud.mig <- c(457:517)
hs.mig <- c(518:578)
ls.mig <- c(579:639)
fam.mig <- c(640:700)

num.agents <-  1166 
aspire.enough <- 507 # aspiring migrants

table(fr.k1.data[,p.all.illegal])
table(fr.k2.data[,p.all.illegal])
table(fr.k3.data[,p.all.illegal])

table(st.k1.data[,p.all.illegal])
table(st.k2.data[,p.all.illegal])
table(st.k3.data[,p.all.illegal])

table(fam.k1.data[,p.all.illegal])
table(fam.k2.data[,p.all.illegal])
table(fam.k3.data[,p.all.illegal])

table(ls.k1.data[,p.all.illegal])
table(ls.k2.data[,p.all.illegal])
table(ls.k3.data[,p.all.illegal])

table(hs.k1.data[,p.all.illegal])
table(hs.k2.data[,p.all.illegal])
table(hs.k3.data[,p.all.illegal])

table(open.k1.data[,p.all.illegal])
table(open.k2.data[,p.all.illegal])
table(open.k3.data[,p.all.illegal])

table(no.k1.data[,p.all.illegal])
table(no.k2.data[,p.all.illegal])
table(no.k3.data[,p.all.illegal])

# Use only data with enforcement level = Fig 2
select.policy.level <-  0.3

df.list <- lapply(df.list, function(x) {
  x <- x[which(x[,p.all.illegal] == select.policy.level),]
  x} )

# Reduce all datasets to the same size for subtraction and comparison

df.list <- lapply(df.list, function(x) {
  x <- x[1:700,] 
  x} )

# Add names
names.df.list <-  Cs(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data,
                     no.k2.data,fr.k2.data,fam.k2.data,ls.k2.data,open.k2.data,hs.k2.data,st.k2.data,
                     no.k3.data,fr.k3.data,fam.k3.data,ls.k3.data,open.k3.data,hs.k3.data,st.k3.data)

names(df.list) <- names.df.list

time.point <- 20

for(i in 1:length(df.list))
{
  df.list[[i]]$legal.raw <- (df.list[[i]][,legal.cols[time.point]])
  df.list[[i]]$legal <-  df.list[[i]]$legal.raw / aspire.enough
  df.list[[i]]$fill <- (df.list[[i]][,fi.cols[time.point]])
  df.list[[i]]$semi <- (df.list[[i]][,sl.cols[time.point]])
  df.list[[i]]$all.illegal.raw <- (df.list[[i]]$fill + df.list[[i]]$semi)
  df.list[[i]]$all.illegal <- df.list[[i]]$all.illegal.raw / aspire.enough
  df.list[[i]]$nm <- 1 - (df.list[[i]]$legal + df.list[[i]]$all.illegal)
  df.list[[i]]$all.migs.raw <- df.list[[i]]$legal.raw + df.list[[i]]$all.illegal.raw
}

# Subtract baseline (open)

b.i <- grep("open", names.df.list)

for (i in 1:7)
{
  df.list[[i]]$legal.adj <- df.list[[i]]$legal - df.list[[b.i[1]]]$legal
  df.list[[i]]$nm.adj  <- df.list[[i]]$nm - df.list[[b.i[1]]]$nm
  df.list[[i]]$illegal.adj  <- df.list[[i]]$all.illegal - df.list[[b.i[1]]]$all.illegal
  df.list[[i]]$legal.migs.adj <- (df.list[[i]]$legal.raw - df.list[[b.i[1]]]$legal.raw) / df.list[[i]]$all.migs.raw
  df.list[[i]]$illegal.migs.adj  <- (df.list[[i]]$all.illegal.raw - df.list[[b.i[1]]]$all.illegal.raw) / df.list[[i]]$all.migs.raw
}

for (i in 8:14)
{
  df.list[[i]]$legal.adj <- df.list[[i]]$legal - df.list[[b.i[2]]]$legal
  df.list[[i]]$nm.adj  <- df.list[[i]]$nm - df.list[[b.i[2]]]$nm
  df.list[[i]]$illegal.adj  <- df.list[[i]]$all.illegal - df.list[[b.i[2]]]$all.illegal
  df.list[[i]]$legal.migs.adj <- (df.list[[i]]$legal.raw - df.list[[b.i[2]]]$legal.raw) / df.list[[i]]$all.migs.raw
  df.list[[i]]$illegal.migs.adj  <- (df.list[[i]]$all.illegal.raw - df.list[[b.i[2]]]$all.illegal.raw) / df.list[[i]]$all.migs.raw
}

for (i in 15:length(df.list))
{
  df.list[[i]]$legal.adj <- df.list[[i]]$legal - df.list[[b.i[3]]]$legal
  df.list[[i]]$nm.adj  <- df.list[[i]]$nm - df.list[[b.i[3]]]$nm
  df.list[[i]]$illegal.adj  <- df.list[[i]]$all.illegal - df.list[[b.i[3]]]$all.illegal
  df.list[[i]]$legal.migs.adj <- (df.list[[i]]$legal.raw - df.list[[b.i[3]]]$legal.raw) / df.list[[i]]$all.migs.raw
  df.list[[i]]$illegal.migs.adj  <- (df.list[[i]]$all.illegal.raw - df.list[[b.i[3]]]$all.illegal.raw) / df.list[[i]]$all.migs.raw
}



# Get the average from all row maxes for raw numbers

k1 <- matrix(nrow = 7, ncol = 3)
k2 <-  matrix(nrow = 7, ncol = 3)
k3 <-  matrix(nrow = 7, ncol = 3)


for (i in 1:7)
{
  k1[i,1] <- mean(df.list[[i]]$illegal.adj) 
  k1[i,2] <- mean(df.list[[i]]$nm.adj) 
  k1[i,3] <- mean(df.list[[i]]$legal.adj) 
}

for (i in 1:7)
{
  k2[i,1] <- mean(df.list[[i + 7]]$illegal.adj) 
  k2[i,2] <- mean(df.list[[i + 7]]$nm.adj) 
  k2[i,3] <- mean(df.list[[i + 7]]$legal.adj)
}

for (i in 1:7)
{
  k3[i,1] <- mean(df.list[[i + 14]]$illegal.adj) 
  k3[i,2] <- mean(df.list[[i + 14]]$nm.adj) 
  k3[i,3] <- mean(df.list[[i + 14]]$legal.adj) 
}


k1 <- as.table(k1[-b.i[1],])
k1 <-  k1 * 100
k2 <- as.table(k2[-b.i[1],])
k2 <- k2 * 100
k3 <- as.table(k3[-b.i[1],])
k3 <-  k3 * 100


## ALL KS

all.ks <-  rbind(k1,k2,k3)
v.k <-  matrix(ncol = ncol(all.ks), nrow = nrow(all.ks))

# Reorder

v.k[c(1,2,3),] <- all.ks[c(1,7,13),]
v.k[c(4,5,6),] <- all.ks[c(2,8,14),]
v.k[c(7,8,9),] <- all.ks[c(3,9,15),]
v.k[c(10,11,12),] <- all.ks[c(4,10,16),]
v.k[c(13,14,15),] <- all.ks[c(5,11,17),]
v.k[c(16,17,18),] <- all.ks[c(6,12,18),]

ks.names <- c("K1","K2","Free Movement, K3","K1","K2","Closed, K3","K1","K2",
              "Family, K3","K1","K2","Low-Skilled, K3","K1","K2",
              "High-Skilled, K3","K1","K2","Student, K3")
# all ks

pdf("FigD5",width=7,height=5)
par(pty="s", mar=c(6, 3, 2, 1), lwd = 0.2)
barplot(t(v.k[,c(1,2)]), horiz=T, las=1, col=c("firebrick1","black"),
        border = T,cex.names=0.8, cex.axis = 1.2, font.axis = 1.2,
        cex.lab = 1.2, xlab = "% of Aspiring Migrants", xlim = c(-60,60), axes = T,
        beside = F, names.arg = ks.names)

barplot(t(v.k[,3]), horiz=T, las=1, col="lightblue",
        border = T,cex.names=1.2, cex.axis = 1.2, font.axis = 1.2, axes = F, add = T,
        beside = F)
dev.off()
